Quality check 03, Block group intersections with community districts
Satya Golla and Steve Simon
Created 2022-03-21
This program examines block groups that only partially intersect with community districts.
Load relevant files
# Important abbreviations:
# bg = block group
# bl = block
# cd = community district
# co = county
# nbd=neighborhood
# tr = tract
library(glue)## Warning: package 'glue' was built under R version 4.1.3library(magrittr)## Warning: package 'magrittr' was built under R version 4.1.3library(sf)## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.7 v dplyr 1.0.9
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.0.1 v forcats 0.5.1## Warning: package 'tibble' was built under R version 4.1.3## Warning: package 'dplyr' was built under R version 4.1.3## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::extract() masks magrittr::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x purrr::set_names() masks magrittr::set_names()library(dplyr)
path_name <- "../data/"Contents of bg
load(glue("{path_name}bg.RData"))
glimpse(bg)## Rows: 1,684
## Columns: 4
## $ bg_id <chr> "201039819003", "201030714004", "201030711041", "201030712051~
## $ bg_name <chr> "Block Group 3", "Block Group 4", "Block Group 1", "Block Gro~
## $ bg_area <dbl> 984975.5, 36602869.0, 18313082.5, 26007918.3, 11448660.8, 243~
## $ geometry <POLYGON [°]> POLYGON ((-94.94756 39.3373..., POLYGON ((-95.18756 3~Contents of bl
load(glue("{path_name}bl.RData"))
glimpse(bl)## Rows: 39,903
## Columns: 4
## $ bl_id <chr> "200910518012010", "200910514003006", "200910515002005", "200~
## $ bl_name <chr> "Block 2010", "Block 3006", "Block 2005", "Block 1009", "Bloc~
## $ bl_area <dbl> 38500.823, 51711.677, 104537.047, 73933.762, 95854.528, 19934~
## $ geometry <POLYGON [°]> POLYGON ((-94.64213 38.9774..., POLYGON ((-94.63966 3~Contents of cd
load(glue("{path_name}cd.RData"))
glimpse(cd)## Rows: 59
## Columns: 4
## $ cd_id <dbl> 106, 108, 113, 102, 129, 116, 114, 101, 105, 103, 107, 109, 1~
## $ cd_name <chr> "East Side", "Old Northeast", "Greater Downtown", "Blue Valle~
## $ cd_area <dbl> 281713850, 118237219, 181920811, 216842678, 412671242, 264282~
## $ geometry <POLYGON [°]> POLYGON ((-94.52337 39.0941..., POLYGON ((-94.50777 3~Contents of co
load(glue("{path_name}co.RData"))
glimpse(co)## Rows: 220
## Columns: 19
## $ STATEFP <chr> "29", "29", "20", "20", "20", "29", "29", "20", "20", "29", "~
## $ COUNTYFP <chr> "083", "011", "073", "043", "157", "103", "117", "039", "147"~
## $ COUNTYNS <chr> "00758496", "00758460", "00485003", "00484991", "00485042", "~
## $ GEOID <chr> "29083", "29011", "20073", "20043", "20157", "29103", "29117"~
## $ NAME <chr> "Henry", "Barton", "Greenwood", "Doniphan", "Republic", "Knox~
## $ NAMELSAD <chr> "Henry County", "Barton County", "Greenwood County", "Donipha~
## $ LSAD <chr> "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "~
## $ CLASSFP <chr> "H1", "H1", "H1", "H1", "H1", "H1", "H1", "H1", "H1", "H1", "~
## $ MTFCC <chr> "G4020", "G4020", "G4020", "G4020", "G4020", "G4020", "G4020"~
## $ CSAFP <chr> NA, NA, NA, "312", NA, NA, NA, NA, NA, "312", NA, NA, NA, NA,~
## $ CBSAFP <chr> NA, NA, NA, "41140", NA, NA, NA, NA, NA, "47660", NA, "21380"~
## $ METDIVFP <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N~
## $ FUNCSTAT <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "~
## $ ALAND <dbl> 1805054746, 1533351028, 2961143636, 1019105691, 1857998463, 1~
## $ AWATER <dbl> 91657543, 12152201, 24145021, 12424174, 7951912, 7307094, 161~
## $ INTPTLAT <chr> "+38.3864909", "+37.5007989", "+37.8793472", "+39.7885021", "~
## $ INTPTLON <chr> "-093.7926278", "-094.3440893", "-096.2417321", "-095.1472253~
## $ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-93.51565 3..., MULTIPOLYGON (((~
## $ label <chr> "Henry", "Barton", "Greenwood", "Doniphan", "Republic", "Knox~Contents of cd-intersections
load(glue("{path_name}cd-intersections.RData"))
glimpse(bg_cd_intersection)## Rows: 1,236
## Columns: 3
## $ bg_id <chr> "290950161001", "290950166002", "290950060001", "2909500340~
## $ cd_id <dbl> 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106,~
## $ bg_prop_in <dbl> 0.8919314491, 1.0000000000, 1.0000000000, 1.0000000000, 0.3~glimpse(bl_cd_intersection)## Rows: 15,748
## Columns: 5
## $ tr_id <chr> "29095002300", "29095006300", "29095002200", "29095002200",~
## $ bg_id <chr> "290950023002", "290950063002", "290950022003", "2909500220~
## $ bl_id <chr> "290950023002014", "290950063002017", "290950022003001", "2~
## $ cd_id <dbl> 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106,~
## $ bl_prop_in <dbl> 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 0.997, 0.001, 0.0~Contents of cd-weights
load(glue("{path_name}cd-weights.RData"))
glimpse(bg_counts)## Rows: 1,236
## Columns: 6
## Groups: bg_id [776]
## $ bg_id <chr> "200910500001", "200910500002", "200910500003", "20091050100~
## $ cd_id <dbl> 209, 209, 209, 209, 211, 209, 211, 211, 112, 112, 112, 111, ~
## $ people_in <dbl> 1, 2, 0, 0, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, ~
## $ units_in <dbl> 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ~
## $ people <dbl> 744, 695, 751, 1313, 1313, 1196, 643, 630, 778, 961, 1210, 9~
## $ units <dbl> 341, 311, 377, 648, 648, 549, 315, 444, 287, 335, 461, 430, ~glimpse(bg_list)## chr [1:776] "290950161001" "290950166002" "290950060001" "290950034004" ...glimpse(bl_counts)## Rows: 41,220
## Columns: 6
## $ bg_id <chr> "290950023002", "290950063002", "290950022003", "29095002200~
## $ bl_id <chr> "290950023002014", "290950063002017", "290950022003001", "29~
## $ cd_id <dbl> 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, ~
## $ prop_in <dbl> 1.000000000, 1.000000000, 1.000000000, 1.000000000, 1.000000~
## $ people_in <dbl> 62, 0, 22, 26, 42, 49, 43, 0, 0, 20, 5, 56, 58, 48, 47, 13, ~
## $ units_in <dbl> 18, 0, 9, 18, 22, 22, 14, 0, 0, 10, 4, 12, 19, 24, 28, 4, 48~glimpse(bl_list)## Rows: 17,547
## Columns: 2
## $ bg_id <chr> "200910533013", "200910533024", "200910533024", "200910534291", ~
## $ bl_id <chr> "200910533013008", "200910533024013", "200910533024008", "200910~Subset co to seven counties
clist <- c(
"20091",
"20103",
"20209",
"29037",
"29047",
"29095",
"29165")
co %>%
filter(GEOID %in% clist) -> co
glimpse(co)## Rows: 7
## Columns: 19
## $ STATEFP <chr> "29", "29", "20", "20", "29", "29", "20"
## $ COUNTYFP <chr> "037", "165", "091", "209", "095", "047", "103"
## $ COUNTYNS <chr> "00758473", "00758537", "00485010", "00485065", "00758502", "~
## $ GEOID <chr> "29037", "29165", "20091", "20209", "29095", "29047", "20103"
## $ NAME <chr> "Cass", "Platte", "Johnson", "Wyandotte", "Jackson", "Clay", ~
## $ NAMELSAD <chr> "Cass County", "Platte County", "Johnson County", "Wyandotte ~
## $ LSAD <chr> "06", "06", "06", "06", "06", "06", "06"
## $ CLASSFP <chr> "H1", "H1", "H1", "H6", "H1", "H1", "H1"
## $ MTFCC <chr> "G4020", "G4020", "G4020", "G4020", "G4020", "G4020", "G4020"
## $ CSAFP <chr> "312", "312", "312", "312", "312", "312", "312"
## $ CBSAFP <chr> "28140", "28140", "28140", "28140", "28140", "28140", "28140"
## $ METDIVFP <chr> NA, NA, NA, NA, NA, NA, NA
## $ FUNCSTAT <chr> "A", "A", "A", "C", "A", "A", "A"
## $ ALAND <dbl> 1804223313, 1087283058, 1226679688, 392739381, 1565698757, 10~
## $ AWATER <dbl> 14963188, 16949451, 16319024, 11801597, 30621016, 28508804, 1~
## $ INTPTLAT <chr> "+38.6464737", "+39.3786900", "+38.8839065", "+39.1153842", "~
## $ INTPTLON <chr> "-094.3545467", "-094.7614765", "-094.8223295", "-094.7630866~
## $ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-94.11966 3..., MULTIPOLYGON (((-94.65596 3..~
## $ label <chr> "Cass", "Platte", "Johnson", "Wyandotte", "Jackson", "Clay", ~Plot functions
# These functions plot various census,
# community, or neighborhood geographies.
add_red_border <- function(mapx, geo) {
mapx +
geom_sf(
data=geo,
aes(),
fill=NA,
color="darkred")
}
add_green_interior <- function(mapx, geo) {
mapx +
geom_sf(
data=geo,
aes(),
fill="lightgreen",
color=NA)
}Get id values for all community districts
cd %>%
tibble %>%
distinct(cd_id) %>%
arrange(cd_id) %>%
pull(cd_id) -> cd_listDraw an overview map
co_points <- st_point_on_surface(co)## Warning in st_point_on_surface.sf(co): st_point_on_surface assumes attributes
## are constant over geometries of x## Warning in st_point_on_surface.sfc(st_geometry(x)): st_point_on_surface may not
## give correct results for longitude/latitude dataco_coords <- as.data.frame(st_coordinates(co_points))
co_coords$NAME <- co$NAME
ggplot()+
geom_sf(
data=co,
color="darkred",
fill=NA) +
geom_text(
data=co_coords,
aes(X, Y, label=NAME)) +
xlab("Longitude") +
ylab("Latitude") -> overview_map
plot(overview_map)
High level view
for (i_cd in cd_list) {
# Draw community district in green atop
# counties outlined in red
cd_subset <- filter(cd, cd_id==i_cd)
ti <- cd_subset$cd_name
overview_map +
ggtitle(ti) +
geom_sf(
data=cd_subset,
aes(),
fill="lightgreen",
color=NA) -> map1
plot(map1)
# draw community districts in green with
# block groups outlined in red.
bg_cd_intersection %>%
filter(cd_id==i_cd) %>%
filter(bg_prop_in > 0.10) %>%
pull(bg_id) -> bg_list
bg %>%
filter(bg_id %in% bg_list) -> bg_subset
ggplot(cd_subset) +
ggtitle(ti) +
xlab("Longitude") +
ylab("Latitude") +
geom_sf(
data=cd_subset,
aes(),
fill="lightgreen",
color=NA) +
geom_sf(
data=bg_subset,
aes(),
fill=NA,
color="darkred") -> map2
plot(map2)
# draw block groups and blocks outlined
# in red atop community district in green.
for (i_bg in bg_list) {
bg_cd_intersection %>%
filter(cd_id==i_cd) %>%
filter(bg_id==i_bg) %>%
pull(bg_prop_in) -> p
if (p < 0.1) next
if (p > 0.9) next
bl_subset <- filter(bl, str_sub(bl_id, 1, 12)==i_bg)
ggplot(cd_subset) +
ggtitle(ti) +
xlab("Longitude") +
ylab("Latitude") +
geom_sf(
data=cd_subset,
aes(),
fill="lightgreen",
color=NA) +
geom_sf(
data=bl_subset,
aes(),
fill=NA,
color="darkred") -> map3
plot(map3)
}
}
for (i_cd in cd_list[5:7]) {
next
bg_n1 <- dim(bg_subset)[1]
bg_n2 <- sum(bg_subset$bg_prop_in==1)
bg_n3 <- sum(bg_subset$bg_prop_in > hi & bg_subset$bg_prop_in != 1)
bg_n4 <- sum(bg_subset$bg_prop_in > lo & bg_subset$bg_prop_in <= hi)
bg_n5 <- sum(bg_subset$bg_prop_in <= lo)
bg_message <- glue(
"{bg_n1} block groups touch or intersect community district\n",
"{bg_n2} block groups lie entirely inside community district\n",
"{bg_n3} block groups lie mostly inside community district\n",
"{bg_n4} block groups lie partially inside community district\n",
"{bg_n5} block groups lie mostly outside community district\n")
cat(bg_message)
plot_cd(cd, i_cd, bg, bg_subset$bg_id, "block groups touching or intersecting")
plot_cd(cd, i_cd, bg, bg_subset$bg_id, "block groups touching or intersecting")
for (i_bg in bg_subset$bg_id) {
bg_cd_intersection %>%
filter(cd_id==i_cd) %>%
filter(bg_id==i_bg) %>%
pull(bg_prop_in) -> bg_prop
if (bg_prop < lo) next
if (bg_prop > hi) next
plot_cd(cd, i_cd, bg, i_bg)
plot_bg(cd, i_cd, bg, i_bg, bl)
bl_cd_intersection %>%
filter(cd_id==i_cd) %>%
filter(bg_id==i_bg) %>%
select(bl_id, bl_prop_in) -> bl_subset
bl_n1 <- dim(filter(bl_cd_intersection, bg_id==i_bg))[1]
bl_n2 <- sum(bg_subset$bg_prop_in==1)
bl_n3 <- sum(bg_subset$bg_prop_in > hi & bg_subset$bg_prop_in != 1)
bl_n4 <- sum(bg_subset$bg_prop_in > lo & bg_subset$bg_prop_in <= hi)
bl_n5 <- sum(bg_subset$bg_prop_in <= lo)
bl_n6 <- bl_n1 - bl_n2 - bl_n3 - bl_n4 - bl_n5
bl_message <- glue(
"{bl_n1} blocks in the block group\n",
"{bl_n2} blocks lie entirely inside community district\n",
"{bl_n3} blocks lie mostly inside community district\n",
"{bl_n4} blocks lie partially inside community district\n",
"{bl_n5} blocks lie mostly outside community district\n",
"{bl_n6} blocks lie entirely outside community district\n")
cat(bl_message)
for (i_bl in bl_subset$bl_id) {
bl_cd_intersection %>%
filter(cd_id==i_cd) %>%
filter(bg_id==i_bg) %>%
filter(bl_id==i_bl) %>%
pull(bl_prop_in) -> bl_prop
if (bl_prop > hi) next
if (bl_prop < lo) next
plot_bl(cd, i_cd, bg, i_bg, bl, i_bl)
}
bg_counts %>%
filter(cd_id==i_cd) %>%
filter(bg_id==i_bg) %>%
data.frame %>%
print
bl_counts %>%
filter(cd_id==i_cd) %>%
filter(bg_id==i_bg) %>%
select(-bg_id) %>%
data.frame %>%
print
}
}